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1 Introduction 



The increasing experimental precision makes it mandatory to compute higher order cor- 
rections to physical quantities. This essentially requires the evaluation of multi-loop in- 
tegrals. However, even at two loops it is very often not possible to solve the integrals 
exactly and one has to rely on approximations. A very powerful approach to get, never- 
theless, reasonable results is based on expansions in small quantities. Recently a method 
was developed for three-loop polarization functions combining expansions from different 



kinematical regions with the help of conformal mapping and Pade approximation |I|J 
An essential ingredient for this procedure are massive integrals with vanishing external 
momentum. For their computation usually recurrence relations are used which are based 
on the integration-by-parts technique @,|[]. These relations connect Feynman integrals 
with different powers of their denominators. In many cases they provide a possibility to 
express an integral with given degrees of the denominators as a linear combination of a 
few so-called master integrals with prefactors which are rational functions of the space- 
time dimension D. The construction of such a procedure is a nontrivial problem even at 
two- loop level ||. At three loops up to now only the case of vacuum integrals with one 
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non-zero mass and various numbers of massless lines has been considered For the 

problems of practical interest the direct application of these equations usually leads to 
intermediate expressions which need several hundred megabyte up to a few gigabyte of 
disk space. 

For the two-loop massless master integral some years ago the recurrence procedure has 
been solved explicitly J7J. The solution, expressed in terms of multiple sums, is also used 
in practice || - - at least for those cases where the exponents of the propagators are not 
too small. Concerning the general multi-loop case in a new approach to implement 
recurrence relations || was suggested. There the factors in front of the master integrals, 
which in the following are called coefficient functions, are considered as independent so- 
lutions of the recurrence relations. In || integral representations for these solutions were 
obtained. It appeared, that in some cases these representations can be expressed in terms 
of Pochhammer symbols. As an example the vacuum integrals with four equal masses and 
two massless lines has been considered, and the efficiency of this approach was demon- 
strated by the calculation of the three-loop QED vacuum polarization. 

In this paper we describe the algorithm, suggested in |J , in more details and discuss some 
peculiarities of its implementation in REDUCE [IIJ and FORM | IT| . In Section |2| the 



general L-loop case is considered. In Section y the derived formulas are specified to three 
loops and explicit solutions are given. The implementation is discussed in Section [| and 
finally Section |5| contains our conclusions. 



2 Basic formulas 



Let us in a first step derive the recurrence relations. The combinatoric structure of these 
relations becomes more transparent if we start with the general multi-loop case, keeping 
in mind, however, the application at three-loop level. Therefore let us consider L-loop 
vacuum integrals with N = L(L + l)/2 denominators (This number of denominators pro- 
vides the possibility to express any scalar product of loop momenta as linear combination 
of the denominators; the diagrams of practical interest which usually have less number of 
denominators, can be considered as special cases with some exponents equal to zero.): 



B (a , D) . %, .... n„, D) = [m ~ m ? I'" I (1> 

where Pi (i = 1, . . . , L) are loop momenta and D a = A l Jpi ■ pj — fi a m 2 (a=l,. . . ,N; here 
and below the sum over repeated indices is understood). For convenience we set m = 1 
in the following. The recurrence relations are obtained by acting with (d/dpi) ■ pk on the 
integrand M: 
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DSlB(n, D) = 2A a kl (I- + /i )Af I d+ B(n, D), (2) 



where I c B{. . . , n c , . . .) = B{. . . , n c - 1, . . .) and I+B(. . ._, n c , . . .) = n c B(. . . , n c + 1, . . .) 
(here, no sum over n c may be performed). The factors A H arise from the fact that the 
scalar products pi ■ pj in the numerator which appear as a result of the differentiation are 
expressed in terms of the denominators D a via: 



p k -p l = A a kl (D a + fi a ) . 
Using the identities 

the recurrence relations (0) can be represented as 

(D — L — l)8iB(n, D) = 2^1^(1" + ^ a )B(n, D), (3) 

where we have exploited that the matrices A a and A a can be chosen to be symmetrical. 
If we denote 

A a kl (l- + ^) = A„, A^I d+ = d il , 
Eq. (^|) will read: 

d il ■ A kl - ^^S k ] B(n, D) = 0. (4) 

Let us now diagonalize these relations with respect to the operators d l . Therefore we 
multiply Eq. (|4]) with det(A), which is the cofactor of the matrix element A k i, and 
sum afterwards over k. Finally we arrive at: 



d il ■ det(A) - 2=±=±(d il det(A)) B(n, D) — 0. 



Our goal is to find the coefficient functions f k (n, D) which relate the given integral -B(n, D) 
with the basic master integrals B(n k , D): 

B(n,D) = f k (n,D)B(n k ,D), 
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where the coefficient functions have to fulfill the initial conditions f l (n k ,D) = 5 l k . The 
index k is used to label different sets of indices n which are fixed by the choice of the 
master integrals. Assuming that the master integrals are algebraically independent one 
can conclude that the functions f k (n, D) should be independent solutions of the recurrence 
relations. This means that as soon as we find a set of independent solutions it is possible 
to construct the desired coefficient functions as linear combinations respecting the given 
initial conditions. To this end it turns out that it is useful to construct an auxiliary 
integral representation for these functions where the operators can be written in the 
forrn^]: I d+ — > d/dxj, 1^ — > x&. Then the differential equation corresponding to Eq. (^) has 
the solution g(x a ) = det(A)^ D_L_1 ^ 2 = P(x a + /i a ) ( ' D ~ L ~ 1 ^ 2 ) where P(x a ) is a polynomial 
in x a of degree L: 



P(x a ) = det(A a kl x a ). 

Therefore it is tempting to consider the "Laurent" coefficients of the function g(x a ) in 
order to determine the coefficient functions: 



f(R, D) = ^ § ■ ■ ■ § |eg£ det(AUx a + ^)) {D ~ L - 1)/2 . (5) 

Here the integral symbols denote N subsequent complex integrations. The contours de- 
pend on the index k and will be specified below. Acting with Eq. (|3|) on (|) leads (up 
to surface terms) to the corresponding differential operator acting on g(x a ), which gives 
zero. The surface terms can be removed by a proper choice of the complex contours: ei- 
ther closed contours or integration paths which end in infinity. For the last case one has 
to consider analytical continuations of D from large negative values. Note that Eq. (|5]) 
is a solution of relation (|2|) and thus the different choices of the contours correspond to 
different solutions which are enumerated with the index k. 

The solutions (||) satisfy by construction the following condition: 



f k (n, D) = P(I- + /i a )/ fc (n, D - 2) . (6) 

In general, however, we are interested in coefficient functions, f k (n, D), which correspond 
to a specific set of master integrals and hence a linear combination of f k (n,D) with 
coefficients depending on D has to be considered. Then Eq. (|6|) gets more complicated as 
a mixing among the different functions is possible: 



f k (n, D) = Sf(D)Pfc + n a )f(n, D - 2), 

2 Note that before the transition to this representation is performed the order of the operators 
has to be reversed. 
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where Sf{D) are rational functions of the space-time dimension D only. 



3 Three-loop case 



In this section we specify the general L-loop formulas derived in the previous one to the 
three-loop case with four equal masses and two massless lines, i.e., \i\ = /i 2 = 0,/i 3 = 
/i4 = /i5 = /U 6 = 1. Then Eq. ([I]) becomes: 



B(n,D) = -. 



2sf"i-3Q /' /' /' d D pd D kd D l 



n"i n"2 n"3 n" 4 n"s n"6 ' 
u 1 u 2 u 3 u 5 ±J a 



with 



D l = k 2 , D 2 = I 2 , D 3 = (p + k) 2 - m 2 , D i = (p + I) 2 - m 2 
D 5 = (p + k + I) 2 - m 2 , D 6 =p 2 -m 2 . 

The coefficient functions from Eq. (|5|) read: 



f(n, D) = ^ ¥ f^...f $P{x u x 2 , x 3 + 1, . . . , x 6 + l) D ' 2 - 2 , (7) 



where the polynomial P(x%, . . . , Xq) is given by 



P(xi, . . . , Xq) = (Xi + X 2 )(xiX 2 - X3X4 - X 5 X 6 ) + (x 3 + X 4 )(-XiX 2 + X3X4 - X 5 X 6 ) 
+ (x 5 + X 6 )(—X 1 X 2 - X3X4 + X 5 X 6 ) + XiX 3 X 6 + X1X4X5 + x 2 x 3 x 5 
+ X 2 X 4 X 6 . 

In this equation we have omitted an overall factor (—1/4) as it would lead to a trivial 
rescaling which cancels after considering initial conditions. 

Let us now compute the coefficient functions which correspond to the choice of the master 
integrals given in HH where B(n, D) has been written in the form: 



B(n, D) = N(n, D)B(0, 0, 1, 1, 1, 1, D) + M(n, D)B(1, 1, 0, 0, 1, 1, D) 
+ T(n, D)B(0, 0,0,1,1,1,D). 

This leads us to the following normalization conditions: 



N(0, 0, 1, 1, 1, 1, D) = 1, N(l, 1, 0, 0, 1, 1, D) = 0, N(0, 0, 0, 1, 1, 1, D) = 0, (8) 
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M(0,0,1,1,1,1,D) = 0, M(l, 1,0, 0,1,1,5) = 1, M(0, 0,0, 1,1, 1,5) = 0, 
T(0, 0, 1, 1, 1, 1, 5) = 0, 5(1, 1, 0, 0, 1, 1, 5) = 0, 5(0, 0, 0, 1, 1, 1, 5) = 1 . (9) 

In a first step the function N(n, 5) is considered. The last two conditions of Eq. (||) are 
satisfied if the contours for the massive indices are chosen to be small circles around zero, 
i.e., Xi = 0. In practice it is useful to perform a Taylor expansion in x 3 ,X4,X5 and x§ 
where, according to the residuum theorem, only one term leads to a result different from 
zero (Here and in the following the proportional sign ("oc") is used as the 5 dependent 
factor is fixed at the end in accordance with the normalization given in Eq. (§).): 



N(n,D)oc 



~TT\ — 772" 

X l X 2 



(n 3 -l)!...(rc 6 -l) 



T yP(x l ,X 2l X 3 + 1, 



X 6 



+ l)D/2-2 



x 3 ,...,x 6 =0 



The remaining integrals over x\ and x 2 , for which in general (because of the derivatives) 
the exponent of the polynomial P(xi, . . . , x 6 ) is reduced, can easily be expressed in terms 
of Pochhammer symbols (a) n = T(a + n)/T(a): 

I I dxy&xg p/„ ^111 i\D/2-2-c (D/2-l)- c (4-3D/2)„ 1+n2+3 c 

IJ^^H^M.Mj OC 4 (n 1+ n 2+ 3 C ) (2 _ D/2)ni+c(2 _ D/2)n2 + c , 

where T(x) is Eulers T function. 

Let us next turn to the function M(n,D) which can be treated in analogy. The only 
difference is that due to the symmetry 5(1,1,0,0,1,1,/)) = 5(1,1,1,1,0,0,1)) it is 
profitable to consider the sum of the solutions with a Taylor expansion performed in the 
variables (xi, x 2 , x 3 , x±) and (xi, x 2 , x$, Xq): 



M(n, D) = m(nx,n 2 , n 3 , n 4 , n 5 , n 6 , 5) + m(n x , n 2 , n 5 , n 6 , n 3 , n 4 , 5) 
where m(n,D) is given by 



m(n, D) oc 



dx^dxR 

"5 «6 

x 5 x 6 



^ni — 1 r^riA — 1 



nU)C d (L-l)l P ( X ^ X ^ X * + 1, • • • , X 6 + l) D/2 - 2 



xi,...,X4=0 



Again the remaining integrals are easily performed with the result 



ftfj|P(0, 0, 1, 1, x 5 + 1, x 6 + lf/i-i-c 

x 5 x fi 

(~)"5+"6 (3-D) n5+ 2c(3-g) n6+ 2c(4-3ZV2) n5+n6+3e 



'OC 



(2-D/2) c (6-2D) n5+nf . 



The case "T" is more complicated. The conditions (0) show that we should find a solution 
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which is non-zero if one of the massive indices is non-positive. It turns out that the 
following combination of "Taylor" solutions leads to the desired result: 

T(ni, n 2 , n 3 , n 4 , n 5 , n 6 , D) = t{n u n 2 , n 3 , n 4 , n 5 , n 6 , D) + t{n u n 2 , n 4 , n 3 , n 6 , n 5 , D) 

+ t(ni, n 2 , n 5 , n 6 , n 3 , n 4 , D) + t(nx, n 2 , n 6 , n 5 , n 4 , n 3 , £>), 

where 



t(n,D) — if n 4 or n 5 or n 6 < 1. (10) 

Four terms are necessary in order to meet the symmetry properties of the initial integral. 
From Eq. (|7|) the following general representation for the solution with the property given 
in Eq. (|T0|) is obtained: 



J J x l x 2 x 3 



X 



„n 4 -l „ng-l „ng-l 
"4 "5 "6 

(n 4 -l)!(n 5 -l)!(ra6-l)! 



P(Xi,X 2 ,X 3 + 1, 



X4,X5,X6=0 



where we have introduced the function t(n, D) as in general the expression given on the 
r.h.s. of Eq. fllTf ) leads to a mixture of the functions N(n, D) and t(n, D). We will moreover 
see below that t(n, D) obeys simple recurrence relations and in addition the extraction 
of t(n, D) and N(n, D) is quite simple. The resulting integral where X 4 — X5 — Xg = is 
given by: 

t(n 1 ,n 2 ,n 3 ,D)=t(n 1 ,n 2 ,n 3 , 1, 1, 1,D) 

7*b/// l^^l- ^ ^ + 1, 1, 1, 1) D/2 - 2 

ffigigt (ag - X!X 2 x 3 + X!X 2 (x! + x 2 - 4)) D / 2 ' 2 . (12) 
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The function t(n\, n 2 , n 3 , D) obeys the following recurrence relations which can be derived 
by using integration-by-parts in Eq. ([12]): 



I+tim, n 2 , n 3 , D + 2) = (D/2 - 1) f-IJIJ + 2lrl 2 + (I 2 ) 2 - il 2 ] t(m, n 2 , n 3 , D) 



I+t{m,n 2 ,n 3 ,D + 2) = {D/2-l) -I^I 3 + 21^ + (l^) 2 - 41^ t(m, n 2 , n 3 , D) 



l^t{nx,n 2 ,n 3y D + 2) 
t(nx,n 2 ,n 3 ,D + 2) 



(D/2 - 1) [2I3 - lrl 2 J t{m, n 2 , n 3 , D), 
I3-) 2 + irlj (-I3 + ir + I2 - 4)1 t(ni, n 2 , n 3 , D) 



(13) 

), 
(14) 

(15) 
(16) 
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The last equation immediately follows from @. 



For n 3 < 1 the contribution from N(n, D) vanishes as the integration around x 3 = 
gives zero and only t(n, D) survives. In this case with the simple change of variables, 
x 3 = + X1X2/2, the integral fll2|) can be reduced to a sum over integrals which 
furthermore can be expressed in terms of Pochhammer symbols. Keeping in mind the 
normalization condition (H) we get the desired solution for t(n,D): 



t(nx, na, n 3 < 1, D - 2c) = 

(D/2-l)- e (D/2-l/2)- c (2--D)(n 1 +n 3 +2c)(2-P)(„ 2+ n 3 +2 C ) 
( _ ) n 3 + C2 (2n 1+ 2„ 2 +6c+3n3) ( 3 /2-D/2) (ni +n3+c) (3/2-D/2) („ 2 + „ 3+c) 

[-13/2] 

V- (l/2) fc (n 3 ) ( - W3 _ 2fc) (I?/2-l/2-c) fc , _s 

2^ (-n 3 -2fc)!(3/2- J D/2+ni+n3+c) fe (3/2-D/2+n 2 +n3+c) fe - " ' J 
k=0 



For n 3 > 1 one can reduce n 3 to 1 with the help of: 



x 



(n3-l)! ( X 3 ~~ x l x 2%3 + x lX 2 (x 1 + X 2 ~ 4)) D/2 2 



(1* 



In connection with the reduction of n 3 two remarks are in order: First, we should mention 
that the reduction of n 3 > 1 to n 3 = 1 can also be expressed through a finite sum which 
is obvious form Eq. ([18]). Second, also for n 3 < 1 the explicit sum in Eq. ([T7|) can be 
replaced by a recurrence procedure using Eqs. ( |15D and ((T6|) in order to arrive at n 3 = 0. 
Afterwards n\ and n 2 are reduced to zero using Eq. (|TTD with n 3 = 0. 



For the case n 3 > 1 one gets after the reduction to n 3 = 1 a set of functions t(ni,n 2 , 1,D — 
2c) with various ni,ri2 and c. The further reduction is done with the help of the relations: 

t(m, n 2 , 1, D) = (2n ( f_"^ 2) [*(ni - 2, n 2 - 1, 1, D - 2) - ±?(ni - 1, n 2 - 1, 0, D - 2)] , 

(19) 

?(m, n 2 , 1, D) = g^g±|?( ni - 1, n 2 + 1, 1, D) + fc^ffi K n 2 + l,0,D), (20) 

which can be obtained from Eqs. (pT3|) - (pT6|) . Note that t(ni, n 2 , 1, D) = t(n 2 , n\, 1, D). 
Therefore these equations can be used to reduce (ni,n 2 ) to (0,-1), (1,0), (0,0) and 
further combinations where n 3 = 0. Then, using Eqs. ( |13D and (|16D for n%, n 2 = —1, 0, 1, 2, 
respectively, we get: 

t(0, -1, 1, D) = |t(0, 0, 1, D) + W0, 0, 0, D), (21) 
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f(l, 0, 1, D) = §§Ef}?(0, 0, 1, D) — ^%%_ A) t(0, 0, 0, D) . (22) 
Finally the shift in the last argument is fixed with the help of 

1(0, 0, 1, D + 2) = - 3( Jjgg,_ 4) [32(£> - 2)7(0, 0, 1, D) + (11 D - 16)t(0, 0, 0, D)\ , 

(23) 

2 )-c n n n) 

4 c (D/2-l-c) (D/2). 



t(0, 0, 0, D - 2c) = a^t 1 * (P /g7 2 Y 2) - c *( > °» 0. ^) » (24) 



which completes the reduction of t(ni, n 2 , n 3 , _D — 2c) to i(0, 0, 1, D) and t(0, 0, 0, D). The 
result for t(ni, n 2 ,n 3 ,D — 2c) is obtained from: 



t(n u n 2 , 1, D) = Q 1 Q — n 2 , 1, D) - 1(0, 0, 1, D)N(n x , n 2 , 1, 1, 1, 1)] , (25) 

which follows from Eq. @. This means that in practice from the final expression of the 
reduction procedure simply the coefficient of t(0, 0, 0, D) has to be extracted in order to 
get the result for t(n,D). Note that due to Eq. fll7|) the functions t(n\, n 2 , 0, D — 2c) 
in Eqs. ([l9| ) - (|23l) are rational in D, i.e., the recursion for n\,n 2 and c are simple one- 
parameter equations of the type F(n) = a(n)F(n — 1) + b(n) where a(n) and b(n) are 
known rational functions. Thus they can easily be solved and expressed through sums 
of Pochhammer symbols. At this point we refrain from listing the resulting formulas 
explicitly. 



4 Implementation 



In this section a possible implementation of the above formulas is discussed. We will 
mainly speak about the most difficult coefficient function ( "T" ) , the others ( "M" and "N" ) 
constitute simple sub-cases and can be considered in complete analogy. The calculation 
of the "T" coefficient function implies subsequent reduction of: 

1) n A ,n b ,n & to 1, 

2) n z to or 1, 

3) m,n 2 to (0,0), (1,0), (0,-1), 

4) D — 2c to D. 

The steps 1) and 2) demand a Taylor expansion of a polynomial raised to non-integer 
power. This can either be done by direct differentiation (see Eqs. ( |TTD and (|I8D), or 
by computing multiple combinatoric sums which arise after rewriting the action of the 
operators d n /dx™. Moreover, the reduction of each parameter can be done independently 
by one of these methods. The steps 3) and 4) can be performed either by using the 
recursion relations fll9D-([2"^]) or by explicitly evaluating the sums. 
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In practice, the quantities of physical interest can be expressed as linear combination of 
vacuum integrals, the number of which can be very large (up to 10 7 ; see, e.g. 0). It is 
very natural to perform the reductions l)-4) not for individual integrals, but for the whole 
expression as simplifications may take place in intermediate steps of the calculation. Of 
course, the method which gives the best performance strongly depends on the computer 
algebra system used. In the following we compare REDUCE 3.6 and FORM 2.3 both 
running on a DEC- Alpha-station with 266 MHz. 



REDUCE. Our experience shows that the method using differentiation in steps 1) and 
2) works about two times faster, although the calculation of the multiple sums allows to 
save computer memory. 

At step 3) the recurrence procedure is also about two times faster than the use of explicit 
sums. The higher demand concerning the memory is not essential at this point as for step 
1) more memory is needed and thus the limitations are given from there. Step 4) needs 
negligible time and memory. Here in practice the recursion approach is used. 



FORM. The explicit differentiation using FORM requires for the variables to be defined 
as "non-commuting" |iT |. This has the consequence that the sorting of the intermediate 



expressions becomes quite slow and needs in addition more main memory. Therefore in 
practice only the method relying on multiple sums are effective for step 1). At steps 2) and 
3) we can choose between the combinatoric-sum and the recursion approach. It turned 
out that their performance is comparable. At step 4) we prefer to use a one-parameter 
table for the functions t(0, 0, 1, D — 2c). 



The implementation in FORM is done in such a way that an expansion in D — 4 is 
performed for the intermediate expressions. Thereby some tricks were used like the appli- 
cation of the "ACCU" function which collects polynomials in D — 4 or the use of tables 
for the Pochhammer symbols |H |. Note that explicit formulas are used for the realization. 
Therefore one has full control on the powers of intermediate poles in D — 4. This actually 
is a problem in the approach based on recurrence relations only. There the expansion in 
D — 4 has to be performed to higher order which has to be fixed after careful examination 
of the recurrence relations. 



It is interesting to mention that the implementation within REDUCE, which even com- 
putes the full .D-dependence, is about three times faster than the one using FORM. An 
explanation for this unexpected fact is the possibility to use direct differentiation within 
REDUCE. It might also be possible that during the calculation essential cancellations 
among long sums containing Pochhammer symbols take place. REDUCE with its abil- 
ity to work with rationals can exploit this fact, however, FORM cannot. We should also 
mention that the attempt to expand in D — 4 in the "REDUCE" version (using "wtlevel" 
declaration [|TTJ) takes more time but is still faster than FORM. 
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Within FORM we are also in the position to compare the recursion approach [Elffl where 



however, in addition to the recurrence relations also tables were used [121, with the one 



described in this paper. As an example we considered the contribution of the non-planar 
diagram to the photon propagator and performed an expansion up to 0(q 12 ) where q is 
the external momentum. It turned out, that the "improved" recursion works better for 
planar vacuum integrals contributing to this diagram, however, worse for the non-planar 
ones. The joined variant which uses "improved" recursion for the planar diagrams and 
the new method for non-planar ones gives a gain of roughly 8 % as compared to the use 
of only the "improved" recursion. 



5 Conclusions 



The new approach described in this paper allows to solve the recurrence relations for 
three-loop vacuum integrals explicitly and expresses the result in terms of multiple sums 
containing Pochhammer symbols. Nevertheless, the final choice of the algorithm which 
leads to the best performance strongly depends on the peculiarities of the problem and 
also on the symbolic language used. 

For the particular problem considered here — three-loop integrals with four massive and 
two massless lines — we get better results using REDUCE 3.6 than using FORM 2.3 from 
point of view of CPU time and in the presence of sufficient memory resources. Certainly, 
FORM stands beyond the competition as soon as the intermediate expressions become 
very large and enormous swapping to the hard disk is necessary. 
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